//
//  Move.h
//  phi_four
//
//  Created by Snir Gazit on 3/9/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#ifndef Move_h
#define Move_h
#include "../HelpLibs/RandomGen.h"
#include "Lattice.h"
#include "../HelpLibs/QmcParameters.h"
namespace ba = boost::accumulators;
typedef ba::accumulator_set<double, ba::stats<ba::tag::mean> > acc;

typedef RandomVariable<MTGen, UniReal> UniRealVariable;
typedef RandomVariable<MTGen, UniInt> UniIntVariable;
typedef RandomVariable<MTGen, Binomial> BinomialVariable;
typedef RandomVariable<MTGen, Discrete> DiscreteVariable;
typedef UniRealVariable::param_type UniRealVariableParams;
typedef UniIntVariable::param_type UniIntVariableParams;
typedef BinomialVariable::param_type BinomialVariableParams;
typedef DiscreteVariable::param_type DiscreteVariableParams;

double calcProb(double lambda, double beta);
using namespace std;
class Move {
protected:
	MTGen *eng_;
	Lattice *lat_;
	UniRealVariable zeroOne_;
public:
	void init(MTGen * eng, Lattice *lat);
	void makeMove();
	bool accept(double p);
	bool checkMove();
};

// Worm Alg Moves
template<class BondData>
class ShiftMove: public Move {
	UniIntVariable chooseDirection_;
//	BinomialVariable createRemove_;
	BoltzWeightType CreateRatio_, RemoveRatio_;
	double g_;
public:
	void init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
			BoltzWeightType boltzRatio);
	void makeMove();
};
template<class BondData>
class JumpMove: public Move {
	UniIntVariable chooseSite_;
	BoltzWeightType boltzRatio_;
public:
	void init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
			BoltzWeightType boltzRatio);
	void makeMove();
};

template<class BondData>
class SwitchMove: public Move {
	BoltzWeightType boltzRatio_;
	double ProbRemove_;
public:
	void init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
			BoltzWeightType boltzRatio);
	void makeMove();
};

// Shift Move
template<class BondData>
void ShiftMove<BondData>::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseDirection_.init(eng,
			UniIntParams(0, BondData::movesPerDim() * DIM - 1));
//	createRemove_.init(eng, BinomialVariableParams(1, 0.5));
	//openCloseRatio_ = qmcParameters.openCloseRatio.getVal();
	g_ = qmcParameters.g.getVal() * (BondData::gFactor());
	RemoveRatio_[0] = 0.;
	for (int j = 0; j < MAX_BOND; j++) {
		CreateRatio_[j] = boltzRatio[j] * g_;
		if (j != (MAX_BOND - 1))
			RemoveRatio_[j + 1] = 1 / (boltzRatio[j] * g_);
	}
}
template<class BondData>
void ShiftMove<BondData>::makeMove() {
// Get Ira
	Lattice::SiteIter ira = lat_->getIra();
// Choose Direction
	int dir = chooseDirection_.draw();
	BondData &cBondData = (ira->getBondArray<BondData>())[dir];
	int &cTotalBonds = (ira->getTotalBonds<BondData>()), cTotalBondsOther =
			ira->getTotalBondsOther<BondData>();
	int &bondVal = cBondData.getBondVal();
	Lattice::SiteIter nSite = cBondData.neighbourSite_;
	if (cBondData.createRemove_ == CREATE) {
		int totalBonds = nSite->getTotalBonds<BondData>();
		int otherTotalBonds = nSite->getTotalBondsOther<BondData>();
		double p = (CreateRatio_[totalBonds + otherTotalBonds]
				* BondData::addRatio(totalBonds, otherTotalBonds))
				/ double(bondVal + 1);
		if (lat_->isClosed())
			p = p * 2;
		else if (nSite == lat_->getMasha())
			p = p * 0.5;
		if (this->accept(p)) {
			bondVal++;
			assert(bondVal < MAX_BOND);
			(nSite->getTotalBonds<BondData>())++;
			lat_->setIra(cBondData.neighbourSite_);
		}
	} else if (bondVal != 0) {
		int totalBonds = cTotalBonds;
		int otherTotalBonds = cTotalBondsOther;
		double p = (RemoveRatio_[totalBonds + otherTotalBonds]
				* BondData::removeRatio(totalBonds, otherTotalBonds))
				* double(bondVal);
		if (lat_->isClosed())
			p = p * 2;
		else if (nSite == lat_->getMasha())
			p = p * 0.5;
		if (this->accept(p)) {
			bondVal--;
			assert(bondVal < MAX_BOND);
			cTotalBonds--;
			lat_->setIra(cBondData.neighbourSite_);
		}
	}
	assert(this->checkMove());
}
// Jump Move
template<class BondData>
void JumpMove<BondData>::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters_, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
//boltzRatio_ = boltzRatio;
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	boltzRatio_ = boltzRatio;
}

template<class BondData>
void JumpMove<BondData>::makeMove() {
	double p;
	Lattice::SiteIter ira = lat_->getIra();
	int &cTotalBonds = ira->getTotalBonds<BondData>();
	int cTotalBondsOther = ira->getTotalBondsOther<BondData>();
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	int &jTotalBonds = jump->getTotalBonds<BondData>();
	int jTotalBondsOther = jump->getTotalBondsOther<BondData>();
	p = (boltzRatio_[jTotalBonds + jTotalBondsOther]
			* BondData::addRatio(jTotalBonds, jTotalBondsOther))
			/ (boltzRatio_[cTotalBonds - 1 + cTotalBondsOther]
					* BondData::addRatio(cTotalBonds - 1, cTotalBondsOther));
	if (this->accept(p)) {
		cTotalBonds--;
		jTotalBonds++;
		lat_->setIra(jump);
		lat_->setMasha(jump);
	}
	assert(this->checkMove());
}

// Switch Move
template<class BondData>
void SwitchMove<BondData>::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters_, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	boltzRatio_ = boltzRatio;
}
template<class BondData>
void SwitchMove<BondData>::makeMove() {
	Lattice::SiteIter ira = lat_->getIra();
	int &cTotalBonds = ira->getTotalBonds<BondData>();
	int &cTotalBondsOther = ira->getTotalBondsOther<BondData>();
	double p = BondData::addRatioOther(cTotalBonds - 1, cTotalBondsOther)
			* BondData::removeRatio(cTotalBonds, cTotalBondsOther);
	if (this->accept(p)) {
		lat_->setSectorType(BondData::getOtherSectorType());
		cTotalBonds--;
		cTotalBondsOther++;
	}
	assert(this->checkMove());
}
// WormLoop Move
typedef boost::array<double, 2 * 2 * DIM> Weights;

template<class BondData>
struct LoopData {
	BondData * bondData;
//	double totalWeight;
//	int totalBonds;
};

template<class BondData>
class CreateLoopMove: public Move {
// typedef typename LoopType::BondData BondData;
	typedef vector<LoopData<BondData> > LoopDataArray;
	typedef typename LoopDataArray::iterator LoopDataArrayIter;
	UniIntVariable chooseSite_;
	BoltzWeightType CreateRatio_, RemoveRatio_;
	Weights weights_;
	LoopData<BondData> currentLoopData_;
	LoopDataArray currentLoop_;
	int avWormLength_;
	double g_;
	acc eraseWorm_;
	int latSize_;
public:
	void init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
			BoltzWeightType boltzRatio);
	void setCoupeling(BoltzWeightType boltzRatio, double g);
	int makeMove();
	void nextStep();
	void eraseLoop();
	void getWeights();
	double getEraseProb();
	int chooseDirection();
	bool continueLoop();
	void setAvWormLength(int avWormLength);
};

template<class BondData>
void CreateLoopMove<BondData>::setCoupeling(BoltzWeightType boltzRatio,
		double g) {
	g_ = g * (BondData::gFactor());
	RemoveRatio_[0] = 0.;
	for (int j = 0; j < MAX_BOND; j++) {
		CreateRatio_[j] = boltzRatio[j] * g_;
		if (j != (MAX_BOND - 1))
			RemoveRatio_[j + 1] = 1 / (boltzRatio[j] * g_);
	}
}

template<class BondData>
void CreateLoopMove<BondData>::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	setCoupeling(boltzRatio, qmcParameters.g.getVal());
	latSize_ = lat_->getNumOfSites();
}

// Create Loop Move 
template<class BondData>
int CreateLoopMove<BondData>::makeMove() {
	int wormLength = 0;
	currentLoop_.clear();
//currentLoop_.resize(avWormLength_);
// Choose First Site
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	lat_->setIra(jump);
	lat_->setMasha(jump);
	getWeights();
	double normFirst = std::accumulate(weights_.begin(), weights_.end(), 0.);
// Start Looping
	do {
		wormLength++;
		nextStep();
		getWeights();
		if (wormLength > 100*latSize_) {
			eraseLoop();
			eraseWorm_(1);
			LoopDataArray().swap(currentLoop_);
			return 0;
		}
	} while (!lat_->isClosed());
	double normLast = std::accumulate(weights_.begin(), weights_.end(), 0.);
	double pc = normFirst / normLast;
	if (!accept(pc)) {
		eraseLoop();
		eraseWorm_(1);
		//	assert(this->checkMove());
		return 0;
	}
//	assert(this->checkMove());
	eraseWorm_(0);
	return wormLength;
}

template<typename BondData>
void CreateLoopMove<BondData>::eraseLoop() {
	Lattice::SiteIter cIra = lat_->getIra();
	for (auto itl = currentLoop_.begin(); itl != currentLoop_.end(); itl++) {
		BondData &cBondData = *(itl->bondData);
		Lattice::SiteIter nextIra = cBondData.neighbourSite_;
		// int &totalBonds = nextIra->totalBonds_;
		int &bondVal = cBondData.getBondVal();
		if (cBondData.createRemove_ == CREATE) {
			bondVal--;
			(nextIra->getTotalBonds<BondData>())--;
		} else {
			bondVal++;
			(cIra->getTotalBonds<BondData>())++;
		}
		cIra = nextIra;
	}
}

template<typename BondData>
void CreateLoopMove<BondData>::getWeights() {
	auto itW = weights_.begin();
// Choose First Site
	auto ira = lat_->getIra();
	const vector<BondData> &cBonds = ira->getBondArray<BondData>();
	int cTotalBonds = (ira->getTotalBonds<BondData>()), cOtherTotalBonds =
			ira->getTotalBondsOther<BondData>();
//	currentLoopData_.totalBonds = cTotalBonds;
	double weightRemove;
	if (cTotalBonds == 0)
		weightRemove = 0;
	else
		weightRemove = (RemoveRatio_[cTotalBonds + cOtherTotalBonds]
				* BondData::removeRatio(cTotalBonds, cOtherTotalBonds));
	for (auto itb = cBonds.begin(); itb != cBonds.end(); itb++) {
		BondData bond = *itb;
		int bondVal = bond.getBondVal();
		if (bond.createRemove_ == CREATE) {
			Lattice::SiteIter nSite = bond.neighbourSite_;
			int totalBonds = nSite->getTotalBonds<BondData>();
			int otherTotalBonds = nSite->getTotalBondsOther<BondData>();
			double weight = (CreateRatio_[totalBonds + otherTotalBonds]
					* BondData::addRatio(totalBonds, otherTotalBonds))
					/ double(bondVal + 1);
			*itW++ = min(1., weight);
		} else {
			double weight = weightRemove * double(bondVal);
			*itW++ = min(1., weight);
		}
	}
}
template<typename BondData>
int CreateLoopMove<BondData>::chooseDirection() {
	double totalWeight = accumulate(weights_.begin(), weights_.end(), 0.);
	//currentLoopData_.totalWeight = totalWeight;
	double draw = totalWeight * zeroOne_.draw(), sum = 0.;
	int chosen = -1;
	for (auto itW = weights_.begin(); sum < draw; itW++) {
		++chosen;
		sum += *itW;
	};
	if (chosen == -1)
		chosen = 0;
	assert(!(chosen < 0));
	assert(!(chosen >= weights_.size()));
	return chosen;
}
template<typename BondData>
void CreateLoopMove<BondData>::nextStep() {
//	chooseDirection_.setParams(DiscreteVariableParams(weights_));
	int chosenDir = chooseDirection();
	Lattice::SiteIter cIra = lat_->getIra();
//    int spaceDir=chosenDir/(2*DIM);
	BondData &cBondData = (cIra->getBondArray<BondData>())[chosenDir];
	Lattice::SiteIter nextIra = cBondData.neighbourSite_;
// int &totalBonds = nextIra->totalBonds_;
	int &bondVal = cBondData.getBondVal();
	if (cBondData.createRemove_ == CREATE) {
		bondVal++;
		assert(bondVal < MAX_BOND);
		(nextIra->getTotalBonds<BondData>())++;
	} else {
		bondVal--;
		assert(bondVal >= 0);
		(cIra->getTotalBonds<BondData>())--;
	}
	lat_->setIra(nextIra);
	currentLoopData_.bondData = &cBondData;
	currentLoop_.push_back(currentLoopData_);
//	assert(this->checkMove());
}
template<typename BondData>
double CreateLoopMove<BondData>::getEraseProb() {
	return boost::accumulators::mean(eraseWorm_);
}
template<typename BondData>
void CreateLoopMove<BondData>::setAvWormLength(int avWormLength) {
	avWormLength_ = avWormLength;
}
#endif
